Oscillatory Climate Modes in the Indian Monsoon, North Atlantic, and Tropical Paciﬁc

Thispaperexploresthethree-wayinteractionsbetweentheIndianmonsoon,theNorthAtlantic,andthetropical Paciﬁc. Four climate records were analyzed: the monsoon rainfall in two Indian regions, the Southern Oscillation index for the tropical Paciﬁc, and the NAO index for the North Atlantic. The individual records exhibit highly signiﬁcant oscillatory modes with spectral peaks at 7–8yr and in the quasi-biennial and quasi-quadrennial bands. The interactions between the three regions were investigated in the light of the synchronization theory of chaotic oscillators. The theory was applied here by combining multichannel singular-spectrum analysis (M-SSA) with a recently introduced varimax rotation of the M-SSA eigenvectors. A key result is that the 7–8-yr and 2.7-yr oscillatory modes in all three regions are synchronized, at least in part. Theenergy-ratioanalysis,aswellastime-lagresults,suggeststhattheNAOplaysaleadingroleinthe7–8-yrmode. ItwasfoundtherewiththattheSouthAsianmonsoonisnotslavedtoforcingfromtheequatorialPaciﬁc,althoughit doesinteractstronglywithit.Thetime-laganalysispinpointedthistobethecaseinparticularforthequasi-biennial oscillatory modes. Overall,theseresultsconﬁrmthattheapproachofsynchronizedoscillators,combinedwithvarimax-rotated M-SSA, is a powerful tool in studying teleconnections between regional climate modes and that it helps identify the mechanisms that operate in various frequency bands. This approach should be readily applicable to ocean modes of variability and to the problems of air–sea interaction as well.


Introduction
Connections between the Southern Oscillation in the tropical Pacific and the Indian monsoon have been explored since the pioneering work of Sir Gilbert Walker and associates in the early twentieth century. More recently, some attention has also been given to the possible effects of the North Atlantic Oscillation (NAO) on the Indian monsoon.
In this paper, we concentrate on potential teleconnections among all three regions: the Indian monsoon, the tropical Pacific, and the North Atlantic. Most teleconnection studies so far were based on correlations between climatic time series from the regions of interest (Wallace and Gutzler 1981); often the correlations found in this way were no larger in absolute magnitude than 0.6, and they did not clarify in which frequency band the teleconnection under investigation was most active.
To explore in greater depth the commonalities between the well-documented records of the NAO index, the Southern Oscillation index (SOI), and the Indian monsoon rain records-which constitute the focus of our paper-we adopt here the viewpoint of synchronization between chaotic oscillators. From this point of view, each regional climatic index may be considered as a chaotic oscillator that represents the regional climate dynamics; teleconnections then arise as a coupling in space between such climatically active regions (Duane 1997;Duane and Tribbia 2004).
Our synchronization analysis is based on singularspectrum analysis (SSA) and it is partly described in Feliks et al. (2010). It differs from that of previous studies that relied on SSA and on multichannel SSA (M-SSA), as reviewed by Ghil et al. (2002b), inasmuch as we use here M-SSA to examine cross-spectral properties of the climate records from different regions. Furthermore, we rely-and improve upon-the recently proposed varimax rotation of M-SSA eigenvectors (Groth and Ghil 2011) in order to enhance the separation of oscillations and the identification of synchronized oscillators. In this context we show here, once more, that, in the absence of rotation, M-SSA is subject to a degeneracy problem and can generate spurious correlations. Feliks et al. (2010) applied some of the present ideas and methods to study teleconnections between the North Atlantic, the eastern Mediterranean, and the Ethiopian Plateau. In all the climatic indices at their disposal, they found prominent oscillatory activity in the 7-8-yr frequency band, as well as synchronization with the NAO index they used. An analysis based on the energy ratios of the oscillatory modes they had identified raised the possibility that the origin of this 7-8-yr mode may lie in the North Atlantic. These authors suggested, therefore, that this mode may arise from changes in the position and the intensity of the Gulf Stream Front (cf. Jiang et al. 1995b;Dijkstra and Ghil 2005;Feliks et al. 2011;Brachet et al. 2012). They also proposed a physical mechanism of moisture flux for the teleconnections between the North Atlantic, Ethiopian Plateau, and the eastern Mediterranean.
The study of Feliks et al. (2010) is extended herein to the much wider areas of the South Asian subcontinent and the Pacific, and it is methodologically improved and sharpened, following Groth and Ghil (2011). The paper is structured as follows. In section 2, we review briefly the concepts and methods of M-SSA and of the varimax rotation of M-SSA eigenvectors; further methodological details are provided in the appendix. The climatic time series we use are presented in section 3, where their univariate spectral properties are also summarized. We study in section 4 the synchronization between the oscillations that were thus identified in section 3. Concluding remarks follow in section 5.

Methodology
Univariate SSA and its multivariate extension to M-SSA rely on the classical Karhunen-Lo eve spectral decomposition of time series. Broomhead and King (1986a,b) introduced them into dynamical system analysis, as a robust version of the Mañ e-Takens idea to reconstruct dynamics from a time-delayed embedding of time series. The method essentially diagonalizes the lag-covariance matrix into eigenvectors and eigenvalues, with the former describing a set of space-time empirical orthogonal functions (ST-EOFs). Projecting the time series onto the ST-EOFs yields the corresponding principal components.
These components provide a ''skeleton'' of the dynamical system's structure; by using linear combinations thereof, we are able to reconstruct trends, oscillatory modes, and irregular ''noise'' (Vautard and Ghil 1989;Ghil and Vautard 1991;Ghil et al. 2002b). The reconstructions associated with the individual ST-EOFs are referred to as reconstructed components (RCs), and each RC represents a narrowband part of the full frequency spectrum. In particular, two RCs associated with nearly equal eigenvalues and dominant frequencies are referred to as an oscillatory pair (Vautard and Ghil 1989;Plaut and Vautard 1994;Ghil et al. 2002b). To avoid misinterpreting random fluctuations as oscillations, a statistical test is performed by comparing the variance captured by such a pair of EOFs with that present for the same pair in a large ensemble of red-noise surrogates (Allen and Smith 1996).
Recently, Groth and Ghil (2011) have demonstrated that a classical M-SSA analysis suffers from a degeneracy problem, with eigenvectors not separating well between distinct oscillations when the corresponding pairs of eigenvalues are similar in size. This problem is a shortcoming of principal component analysis in general (Richman 1986;Jolliffe 2002), not just of M-SSA in particular. To reduce mixture effects and to improve the physical interpretation of the results, Groth and Ghil (2011) have proposed a subsequent varimax rotation of the ST-EOFs. To avoid a loss of spectral properties (Plaut and Vautard 1994), their proposed algorithm has been modified somewhat from the common varimax rotation of Kaiser (1958), by taking into account the spatiotemporal structure of ST-EOFs (see the appendix for further details).

The climatic records and their oscillatory modes
In this section, we study the spectral properties of the individual climatic records by means of a single-channel SSA analysis. The climatic records are provided as monthly averages, and they are dominated by strong seasonal activity. Since we are interested in lowerfrequency activity, on interannual time scales, we first remove all frequencies higher than 0.5 cycles per year and proceed with an annual sampling rate. Such lowpass filtering avoids in particular aliasing effects in the annually subsampled records.
We applied a Chebyshev type I filter that has a very steep rolloff at the desired cutoff frequency and a nearly constant magnitude response in the passing band. When compared to a 12-month (i.e., an annual mean) moving average or to a 3-month one (i.e., a seasonal mean), the Chebyshev filter is much closer to an ideal low-pass filter and does practically not alter the spectral properties in the desired low-frequency band (see the appendix). In the filtering process of the monthly record, we pass it in both the forward and reverse direction in order to ensure zero phase distortion (Oppenheim and Schafer 2009). Figure 1 shows the low-pass filtered climatic records, from which we finally derived annually sampled time series by simply taking all July values. Single-channel SSA is applied to each of these time series, and the statistical significance is assessed by means of Monte Carlo SSA (MC-SSA; Allen and Smith 1996), with an ensemble size of 1000 surrogate time series. To check the robustness of the peaks identified via SSA, the Mann and Lees (1996) version of the multitaper method (MTM) is applied, with three tapers [see also Ghil et al. (2002b) and the documentation available at http://www. atmos.ucla.edu/tcd/ssa]. As for SSA, we assess the statistical significance of the MTM results against a rednoise null hypothesis. A summary of the modes identified at a 95% level of significance-in both MC-SSA and MTM, respectively-are given in Table 1, along with the percentage of variance captured by the corresponding SSA mode.

a. SOI
The SOI is defined as a suitably normalized difference in the month-to-month fluctuations in the mean sea level pressure between Tahiti and Darwin. We use the monthly index between years 1876 and 2006, as compiled by the Australian Bureau of Meteorology (see http://www.bom.gov.au/climate/current/soi2.shtml).
After a reduction of the sampling rate to annual values, as explained in section 2, we applied SSA with a window width of M 5 40 yr; the results in Table 1 were found to be robust with respect to reasonable variations in M. The oscillatory modes that are statistically significant can  be classified into four broadband peaks, with a period of 6-7, 5, 3.6, and 2.5 yr, respectively.
The 2.5-yr band contains 7% of the total variance and corresponds to the well-known quasi-biennial oscillation component of ENSO (Rasmusson and Carpenter 1982;Jiang et al. 1995a;Palu s and Novotn a 2006). The 3.6-yr peak contains 10% of the variance and is probably identifiable with the low-frequency or quasi-quadrennial oscillation of ENSO (Jiang et al. 1995a;Ghil and Robertson 2000).

b. Indian monsoon rainfall
The Indian summer monsoon rainfall occurs during the months of June-September (JJAS), and its statistics are divided into several homogeneous monsoon rainfall regions by the Indian Institute of Tropical Meteorology (Parthasarathy et al. 1995). Azad et al. (2010) suggested a different partition into homogeneous regions, based on spectral properties, but here we follow the standard one. These regions differ in the onset time of the monsoon, in the total amount of rain, and in the mechanisms that influence the rainfall. In this study, we analyze records from the core and the peninsular regions that cover 138 yr, from 1872 to 2008 (see Figs. 1b,c for the low-pass filter version of these two records). The core region represents the central and western parts of India, while the peninsular region represents its southern part. As statistically significant oscillatory modes, we can mainly identify three broadband peaks of 2.2-2.8, 3.2-3.6, and 7.3-7.6 yr. In the peninsular record, we also obtain a significant peak at 5 yr (cf. Table 1).

c. NAO
The NAO has been extensively studied over the last two decades (Hurrell et al. 2003), and its strength has been quantified using several indices. We use here the monthly index over the years 1823-2008, as compiled by Jones et al. (1997) (see Fig. 1 for the monthly low-pass filtered record). This index is the difference between the normalized sea level pressure at Gibraltar and the normalized sea level pressure over southwestern Iceland; Jones and colleagues used early instrumental data back to 1823. The dataset has been modified in November 2000, and the effect of this change is most evident in the summer (see http://www.cru.uea.ac.uk/cru/data/ nao/). The significant oscillatory modes have periods of 7.8, 5.9, 5, 4.3, and 2.7 yr, respectively. The most prominent oscillatory mode, with a period of 7.8 yr, appears in the two leading RCs, RCs 1 and 2, in good agreement with previous findings (Robertson 2001;G amiz-Fortis et al. 2002;Palu s and Novotn a 2004). A quasi-biennial component is present, too, in the closely related Arctic Oscillation, constructed from hemispheric sea level pressures (Trenberth and Paolino 1981;Robertson 2001).

Synchronized modes
We have so far found significant oscillatory modes on interannual time scales to be present in records from the North Atlantic, tropical Pacific, and the Indian subcontinent. Such a coincidence of spectral features in different regions of the world raises the possibility of teleconnections between the climate variations in the regions under consideration.
To study this possibility in detail, we examine next the synchronization of pairs of records by means of M-SSA and look for common oscillatory modes and their role in the dynamics of each region. We analyze synchronization of pairs of records by applying a two-channel M-SSA to the time interval of overlap between the two records of the pair: each channel corresponds to an instrumental or proxy record, nondimensionalized by dividing it by its standard deviation.
For a given oscillatory mode of the M-SSA analysis, we follow Feliks et al. (2010) and evaluate the energy balance between channels by the ratio « of the variances of the two RC components (see the appendix for details). This ratio tends to « 5 1 when the two channels are well balanced and the coupling is strong, and it tends to either « 5 0 or « 5 ' as one process dominates in the oscillatory mode. Note that we rely here on normalized time series, although other physically motivated transformations of the raw indices are possible. Our choice is one of simplicity and, given this choice, the covariance matrix is transformed into a correlation matrix.
One might conjecture that, say when « ( 1 and thus x 2 dominates the mode under study, it is likely to drive the mechanism that activates the mode, because of its much larger energy (Feliks et al. 2010). To buttress further this conjecture, we also study here the time lag t between the RC components to obtain additional information about the leading process in that mode.

a. SOI and NAO
The 130-yr interval of overlap between the SOI and the NAO index is 1876-2006. At first, we carry out an M-SSA analysis with a window width of M 5 40 and without varimax rotation. The resulting estimate of power spectral density (PSD) for the 10 leading ST-EOFs is shown in Fig. 2.
We observe a grouping of the EOFs into pairs with similar PSD characteristics, referred to as oscillatory pairs (see also section 2 and the appendix). These pairs are the data-adaptive equivalent of sine and cosine pairs in Fourier analysis and characterize amplitude-modulated oscillatory modes. We see that their PSD is not always unimodal and can exhibit complex, multimodal characteristics. Still, one typically assigns a dominant frequency to each of these EOFs-for example, derived from the maximum of the PSD (Vautard and Ghil 1989)-in order to simplify the interpretation of the eigenvalue spectrum. Allen and Smith (1997) have discussed the problem of frequency separation in the presence of large colored noise. They proposed an improvement in ascertaining statistical significance of oscillatory modes in singlechannel SSA analysis by estimating a red-noise process from the noisy data. In this paper, however, we do not make any explicit assumption about the nature of the noise in the underlying system and show an improvement in frequency separation for the case of multichannel SSA.
The need for such an improvement is apparent from From the maximum of the PSD, we see that this mode has much less energy in the NAO (heavy solid) than in the SOI (light solid). This low energy of the mode in the NAO index is in agreement with a single-channel SSA of the latter, in which we could not identify such a mode (cf .  Table 1). Hence, this oscillatory mode appears to be significantly present only in the tropical Pacific.
For the remaining EOFs in Fig. 2, the picture is less clear, with overlapping PSDs and multimodal, broadband characteristics of the spectra. Furthermore, when increasing the window width to a value of M 5 50, the isolated 3.6-yr mode in Fig. 2 is also mixed in with other modes (not shown).
A possible approach to overcome this mixing problem is a subsequent M-SSA analysis on a subset of RCs. This means that we exclude EOFs 3 and 4 from the M 5 40 analysis in Fig. 2 and reconstruct the two time series using only RCs 1, 2, and 5-8, for example. Next, we apply M-SSA to this reconstruction but do not constrain the new ST-EOFs to be orthogonal to the previous EOFs. Such an iterative procedure should be able to better separate between different frequencies in a reduced subspace, at least when the variances of two or more pairs are not too close in value.
In general, this mixing problem is due to similar variance in different frequency bands, and it is much more severe when the dimensionality of the system is larger; one needs therewith a more systematic way to disentangle the degeneracy. Groth and Ghil (2011) have proposed to solve this common degeneracy problem in M-SSA by applying a simple-structure rotation (Richman 1986;Jolliffe 2002) of the leading ST-EOFs, while giving due consideration to their spatiotemporal character.
To reduce the risk of an overrotation (O'Lenic and Livezey 1988), we multiply each EOF by its corresponding singular value prior to the modified varimax rotation, as described in greater detail in the appendix. This scaling prevents, in particular, the rotated EOFs from being too localized in space and frequency and the variance spectrum from being too flat. At the same time, the choice of an appropriate number S of rotated EOFs becomes less critical and the resulting components are more consistent over a wide range of S values, as we will see in the following. The price to pay for this scaling is that the rotated EOFs are no longer the eigenvectors of the covariance matrix and a certain variance reduction occurs in the leading EOFs. Still, scaled varimax rotation helps reduce the problem of an artificial variance compression in the eigenvalue spectrum: that is, of a spurious inflation of variance in the largest eigenvalues (Allen and Smith 1996;Groth and Ghil 2011).
The PSD estimate of the leading ST-EOFs after rotation is shown in Fig. 3. The variance localization in frequency is much sharper than prior to rotation, with a clear tendency to unimodal, narrowband PSD features for each ST-EOF. Note that, at the same time, we have slightly increased the window width to a value of M 5 50, in order to better separate nearby oscillations. On the other hand, M values that are too large will tend to split a complex process into an unreasonably high number of narrowband components.
To better understand the frequency separation process through modified varimax rotation, as well as its reliability and reproducibility with respect to parameter changes, we analyze the leading 7.8-yr mode in greater detail. Figure 4a shows the energy ratio « between SOI and NAO of the leading 7.8-yr mode as a function of the relative number of rotated EOFs.
As the number S of rotated EOFs increases, the energy ratio decreases, indicating that the separation in frequency between the two time series increases. However, this separation occurs in three distinct steps, rather than gradually and uniformly: (i) « drops precipitously at The first abrupt drop can be clearly attributed to frequency separation of the leading EOFs. It is especially true for large M values (i.e., M 5 60 and M 5 50) that the EOFs are susceptible to degeneracy prior to any rotation.
In the region (ii) of nearly constant « 5 «(S), the energy ratio is pretty flat, independently of the M value, and it thus gives a robust estimate of the energy balance. In the leading 7.8-yr mode, the SOI has therewith much less energy than the NAO, especially when compared with the results of a classical M-SSA analysis without rotation. This example illustrates quite well the efficacy of varimax rotation in reducing the degeneracy problem and thus the risk of spurious correlations.
The low energy ratio thus obtained suggests that the NAO strongly dominates this oscillatory mode but that its impact in the tropical Pacific Ocean is rather weak and unclear. For large window width (M 5 60), it even turns out that « decreases further beyond S(DM) 21 ' 0.4 and would thus appear to negate any link between SOI and NAO in the 7-8-yr frequency band. This apparent uncoupling could, however, be due to spectral splitting of a broader-band phenomenon.
To confirm our success in solving the degeneracy problem in the 7.8-yr mode, we further analyze the unimodality of the corresponding RCs. We reconstruct separately the two time series by means of the ST-EOF pair and perform a subsequent M-SSA analysis on the resulting RCs, by using the same M value but without rotation. For RCs that are truly unimodal, we expect to find almost the entire energy concentrated in the first two eigenvalues. In the absence of rotation, however, the amount of energy is much lower, and therewith it clearly points to the presence of mode mixing in the leading EOF pair (Fig. 4b). Only rotating a sufficient number of components, S(DM) 21 * 0.1, solves this degeneracy problem and the RCs become unimodal. This result demonstrates quite clearly and convincingly that the frequency separation process is linked indeed to a sudden jump in the energy ratio (see Fig. 4a). Figure 5a shows the complete eigenvalue spectrum, with the statistical significance level. We find three significant oscillatory modes with a period of approximately 7.8, 5.8, and 3.6 yr, respectively; each of these pairs of eigenvalues stands clearly out of the red-noise null hypothesis of MC-SSA (Allen and Smith 1996;Ghil et al. 2002b). Note that, in contrast to the previous univariate SSA analysis, we have slightly increased here the required level of significance to 98%; this increase is meant to compensate for a loss of statistical power of Monte Carlo SSA when the number of channels increases (cf. A. Groth and M. Ghil 2013, unpublished manuscript).
The corresponding RCs are shown in Figs. 5b-d. Note that, although the SSA analysis has been performed on the annual records, the RCs are plotted-here and in Figs. 6-8-using a monthly sampling rate, in order to better visualize high-frequency oscillations, such as the quasi-biennial oscillation. Moreover, the monthly resolution provides a more precise estimation of the time lag t between the two RCs of a pair. To interpolate the annually sampled RCs to monthly resolution, they are first zero padded between the samples and then low-pass filtered with the same Chebyshev filter that has been employed in the annual subsampling (see section 3 and the appendix). The energy ratio, on the other hand, has been directly derived from the annually sampled RCs.
In the leading, 7.8-yr mode of Fig. 5b, we see a clear phase opposition between the SOI and NAO: it is the negative phase of the NAO that is associated with an enhanced SOI index, even though the effect on the latter is quite weak. In the next, 5.8-yr mode (Fig. 5c), the two records are much closer in energy, although this mode is still dominated by the NAO. Even so, it is hard to identify a driving mechanism in this frequency band, given the fact that the approximate 2-yr lead of the NAO with respect to the SOI can easily be changed into a 1-yr lag of the negative NAO.
The 3.6-yr mode (Fig. 5d) is dominated by the SOI, as can be seen from the high value of «. This SOI dominance is also confirmed by single-channel SSA analyses (cf . Table 1), where this mode has only been identified in the SOI and not in the NAO.
The teleconnection between the North Atlantic and the tropical Pacific is of rather moderate strength in our careful, varimax-rotated M-SSA analysis, even though it is supported by mechanisms outlined in previous studies. For instance, Chiang and Vimont (2004) have provided evidence for the extratropical atmosphere influencing the tropical Pacific through a meridional mode, analogous to that present in the tropical Atlantic, where variations in the trade winds in the northern subtropics influence tropical sea surface temperatures (SSTs). Thompson and Lorenz (2004) have furthermore identified relationships between the extratropical annular modes and the circulation in the tropical troposphere; these annular modes, however, were filtered out herein by the annual subsampling (cf. the appendix and Fig. A1).

b. Indian monsoon rainfall and NAO
The years 1871-2008 are the 137-yr interval of overlap between the Indian monsoon and the NAO index. From Table 1 we see that there are fewer possible candidates of common oscillatory modes than in the previous case of SOI and NAO. Hence the degeneracy problem is less pronounced for these two pairs of indices: namely, NAO paired with either the core Indian rainfall or the peninsular rainfall.
First, we analyze the core region and the NAO index. In this case, although the mixing effect is not totally negligible, the significant oscillatory modes in Fig. 6 exhibit practically no multimodal characteristics, even in the absence of rotation (not shown). The significant oscillatory modes in Fig. 6a are in good agreement with the modes that are significant in the individual single-channel SSA analyses of both indices (cf. Table 1). As for the SOI and NAO before, there is again an NAO-dominated 7.7-yr mode and we are led once more to suspect that it originates in the North Atlantic (Fig. 6b). Still, the importance of this mode for the core rainfall is unclear, given its rather weak energy in the latter. Note that, in the single-channel SSA results for the core rainfall, the corresponding mode is only significant at a lower level.
We observe further a quasi-biennial mode in EOFs 3 and 4, with a period of 2.7 yr. In this mode, both records exhibit practically the same energy and no time delay (Fig. 6c), which appears to imply a strong interaction between the two systems in this frequency band.
Next, we carry out the combined M-SSA analysis of the peninsular region and the NAO index. The three significant oscillatory modes in Fig. 7 are in good agreement with the significant modes from the individual, singlechannel SSA analyses summarized in Table 1. Since the variance of the leading 7.6-yr oscillatory mode clearly stands out above the other oscillatory modes, there is no mode mixing and we can exclude EOFs 1 and 2 from the subsequent rotation.
The 7.6-yr mode is shown in Fig. 7b; in it, the energy ratio is close to unity, thus suggesting strong coupling, while the amplitude of the NAO index is only marginally larger than that of the peninsular rain. It follows that the direction of influence is not clear from the energy analysis alone. On the other hand, we do find a lead of the (negative) NAO by 0.5 yr in Fig. 7b that suggests the NAO drives this mode in the Indian peninsular rainfall, which is consistent with the physical arguments below. The amplitude of this joint mode decreases gradually throughout the time interval of record and it vanishes entirely at the end of the interval. This damping is further investigated in each time series separately by single-channel SSA analysis.
In the peninsular rainfall, the leading oscillatory mode is captured by RCs 1 and 2 and it exhibits a period of 7.4 yr, whereas in the NAO index the leading oscillatory mode is still given by RCs 1 and 2 but it has a slightly different period, of 7.7 yr. The reconstructions of the 7.4-yr mode in rainfall (dashed) and of the 7.7-yr mode in the NAO index (solid), respectively, are given in Fig. 8.
From this figure, it follows that neither mode individually is damped since 1950: in fact, the rainfall amplitude is pretty constant throughout the interval, while the amplitude of the NAO index increases from 1950 on. Thus, the damping of the joint mode in Fig. 7b cannot really be attributed to a damping of the two indices in this frequency band. Furthermore, the phase difference in Fig. 8 shifts from phase opposition until 1950 to an inphase behavior after 1985. This phase shift, rather than the amplitudes of the individual modes, results in the damping of the joint mode after 1950 and its almost vanishing amplitude after 1985. To examine the mechanism that might be responsible for NAO effects on the peninsular rainfall, we consider next the spatiotemporal behavior of relevant climate fields. Figure 9 shows global maps of the Pearson productmoment correlation coefficient for the NAO component of the 7.6-yr joint mode in Fig. 7b with each of the following three fields: 200-hPa geopotential heights (Z200; Fig. 9a), mean sea level pressure (MSLP; Fig. 9b), and the SST field (Fig. 9c). The correlations are calculated over the time interval 1876-1940, when the covariability is strongest; each monthly data field was prefiltered by using a Chebyshev type I filter, as done in all previous calculations.
For consistency, we use July values for the three fields, as we do for the NAO index. The atmospheric variables are taken from the twentieth-century reanalysis of Compo et al. (2011), and the SSTs are from Kaplan et al. (1998). For a sample size of 65/2 yr, with the number of years halved to account for the impact of the filtering, a correlation magnitude of 0.35 is statistically significant according to a two-tailed Student's t test.
The MSLP correlation map shows a significant NAO dipolar pattern over the North Atlantic, as expected. It is notable, however, that a wave train arcs southeastward from the North Atlantic at 200 hPa (Fig. 9a), with a significant trough over Central Asia. This trough is centered north of the Tibetan plateau and it coincides with the positive phase of the NAO. Such a trough would tend to weaken the monsoonal anticyclone at upper levels, thus weakening the monsoonal circulation (Saeed et al. 2011). At the surface, MSLP correlations (Fig. 9b) are generally positive over Asia and Africa, although they are weak. The SST anomaly correlations (Fig. 9c) exhibit an El Niño-like pattern in the tropical Pacific, together with NAO-like correlations over the North Atlantic. The former would also tend to be associated with a weaker Indian monsoon. The correlations are not significant if the 1941-2006 period is taken in Fig. 9, suggesting that this physical teleconnection mechanism-from the North Atlantic to Central Asia and on to the Indian subcontinent-seems to be breaking down in the second half of the twentieth century.
After a rotation of the remaining EOFs 3-30, the 2.7and 4.7-yr oscillatory modes in Figs. 7b,c have a clear unimodal behavior too. The energy ratio in the 4.7-yr mode (Fig. 7c) implies that this mode is limited to the peninsular region, having practically no energy in the NAO. In the 2.7-yr mode (Fig. 7b), however, it is the NAO that has more energy, although the evidence for strong coupling in the latter is still weaker than for the 7.6-yr mode.
The teleconnection between the North Atlantic and the Indian summer monsoon (ISM) was studied recently by several authors. Syed et al. (2011) argued that there are two modes in the North Atlantic that influence the ISM. One is the summer circumglobal teleconnection (CGT), as suggested by Ding and Wang (2005); the other one is the summer NAO (Folland et al. 2009). The spatial patterns of these two modes are similar over the North Atlantic.
Different mechanisms have been suggested to explain the role of the CGT and NAO in influencing the summer monsoon in India and Pakistan. Saeed et al. (2011) showed that the positive anomaly at 200 hPa (positive CGT) over western Asia is associated with low surface pressure over northern India and Pakistan. This increases the advection of moist air from the Arabian Sea and the convergence of this marine air inflow over India and Pakistan. Hence the updraft and rainfall increase over Southern Asia for positive CGT, and the opposite occurs for negative anomalies at 200 hPa: that is, negative CGT (cf. Fig. 9a). These authors also claim that the role of the NAO over Central Asia is similar to that of CGT. Ding and Wang (2005) showed that the CGT pattern is most pronounced in August, when the strongest linkage between the Atlantic and Central Asia exist. The pattern in September is very similar to the August one, while in July it is rather different. In the present analysis, we only examined the teleconnection between the North Atlantic and India through the NAO. Rajeevan and Sridhar (2008) also showed that positive SST anomalies over the North Atlantic shift the jet stream northward and the associated circulation changes the upper troposphere's influence on the Indian monsoon rainfall. The 7-8-yr and 2.7-yr oscillation in the troposphere over the North Atlantic may be induced by the same-period oscillation found in the SST field there, as shown by Feliks et al. (2011) andBrachet et al. (2012).

c. SOI and Indian core rainfall
The India core rainfall record and the SOI overlap for 130 yr . The eigenvalue spectrum in Fig. 10a suggests common modes in the quasi-quadrennial and the quasi-biennial frequency band.
The quasi-biennial band is separated into two nearby oscillatory modes. Both the 2.8-yr oscillation in EOFs 1 and 2 and the 2.3-yr oscillation in EOFs 3 and 4 show a lead of the Indian core rainfall with respect to the SOI. This result is consistent with the finding of Rajeevan and Pai (2007) that the South Asian monsoon is not forced by the equatorial Pacific but rather that the two are mutually interacting systems. The 3.6-yr mode in EOFs 5 and 6 ( Fig. 10d) is dominated by the SOI and probably originates in the tropical Pacific.
The interannual oscillations found in the Indian summer rainfall and in the teleconnection with ENSO, as well as with the tropical Pacific SSTs, were the subject of many studies, going back to the early twentieth century (Philander 1990;Sarachik and Cane 2010). Overall, low monsoon rainfall is associated with El Niño yearswarm eastern tropical Pacific and negative SOI-while high monsoon rainfall is associated with La Niña yearscold SSTs off the coast of Peru and positive SOI (Sikka 1980;Rasmusson and Carpenter 1983;Parthasarathy et al. 1988;Kumar et al. 1999;Terray and Dominiak 2005;Rajeevan and Pai 2007)-although there are years during which the opposite is observed. Since the two quasi-biennial oscillations are quite close to each other in an otherwise broadband spectral behavior of the regional climate systems, we examine the possibility of varimax-rotated M-SSA to separate the two. Figure 11a shows the energy ratio between the SOI and the Indian core rainfall in the leading 2.8-yr oscillatory mode, as a function of the number S of rotated components and the window width M.
It turns out that, especially at small M, the EOFs are degenerate in the absence of rotation, with multimodal characteristics of the corresponding RCs (Fig. 11b). As the window width M and thus the spectral resolution increases, this mode mixing vanishes and a more refined solution emerges. Still, this refinement raises immediately the question of what a reasonable spectral resolution might be, given the presence of high stochasticity in the system's dynamics. Nonetheless, the energy ratio is rather constant over a wide plateau of S values and for different M values (Fig. 11a). This robustness is a strong indicator for the reliability of the identified energy ratio of « ' 0.5 and therewith of a genuine separation of the quasi-biennial band into two close but separate modes.

Concluding remarks
In this study we have shown that interannual oscillatory modes exist in the climatic records from India, the North Atlantic, and the tropical Pacific. The records we analyzed were the monsoon rainfall for the Indian core region and the peninsular region for India, the SOI for the tropical Pacific, and the NAO index for the North Atlantic.
Statistically significant oscillatory modes having periods of 2.2-2.8, 3.2-3.6, 5, and 6-8 yr were found in climatic records from two or three regions. Such a clustering of periodicities raises the possibility of teleconnection between these regions. We applied systematically to the study of climatic teleconnections between the regions an approach based on synchronization between chaotic oscillators. This approach was introduced in Feliks et al. (2010) and has been refined here by the technical improvements introduced by Groth and Ghil (2011).
The 7-8-yr and 2.7-yr oscillatory modes in all three regions are synchronized, at least in part. Such an inference is also supported by the recent results of De Viron et al. (2013), who found that four leading modes capture over two-thirds of the variability of 25 climate indices of different physical nature and from different regions, over the last six decades.
Taken together, the energy ratio and the time-lag analysis suggest that the NAO plays a leading role in the 7-8-yr mode. The two oscillations in the North Atlantic-at 7-8 yr and at 2.7 yr-are probably due to oscillations of similar period in the position and strength of the Gulf Stream's SST front. Feliks et al. (2011) showed that statistically significant oscillations in the SST in the 7-10-yr band and at 2.8 yr exist in two regions along the Gulf Stream track: off Cape Hatteras and in the Grand Banks region.
The sharp SST fronts associated with the Gulf Stream spin up an atmospheric jet over the North Atlantic and give rise to oscillations in the 7-10-yr band, similar in period to those observed in the Gulf Stream front itself. In addition, a close correspondence is found between interannual spectral peaks in the observed NAO index and the SST-induced oscillations in an idealized, quasigeostrophic atmospheric model driven by the observed SSTs (Feliks et al. 2011). In particular, significant oscillatory modes with periods of 8.5, 4.2, and 2.8 yr were found by these authors in both the observed and modelsimulated indices.
These modes were shown to be highly synchronized and of similar energy in both observed and modelsimulated time series. This result-together with many previous studies of the 7-8-yr oscillation in SSTs and in sea level pressures over the entire basin (Deser and Blackmon 1993;Moron et al. 1998;Joyce et al. 2000;FIG. 11. MC-SSA analysis of the SOI and the Indian core rainfall, as a function of the relative number S(DM) 21 of rotated EOFs and for different values of the window length M; here D 5 2 is fixed. (a) Energy ratio « of the leading 2.8-yr oscillatory mode and (b) the unimodality index of the corresponding RCs.
Costa and Verdi ere 2002)-shows that the energy flow occurs in the two directions, from the ocean to the atmosphere and back, while Wang et al. (2004) showed that it is the Gulf Stream region alone that affects the NAO.
The physical mechanism of the 7-8-yr variability in the Gulf Stream, in turn, can be credibly linked to an oscillatory gyre mode of the North Atlantic's winddriven circulation (Jiang et al. 1995b;Ghil et al. 2002a;Dijkstra and Ghil 2005;Simonnet et al. 2005;Sushama et al. 2007). The results in the present paper, together with those of Feliks et al. (2010), show that the 7-8-yr oscillatory mode is found around the globe, including over the eastern Mediterranean, East Africa, South Asia, and the tropical Pacific, and that it probably arises in the North Atlantic.
The teleconnection between the North Atlantic and the Indian summer monsoon has been attributed to two separate modes: the summer CGT (Ding and Wang 2005) and the summer NAO (Folland et al. 2009). Different mechanisms have been suggested to explain the role of the CGT and the NAO in affecting the summer monsoon in India and Pakistan. Here we have only examined the role of the NAO, and shown that a strong dipole over the North Atlantic is highly correlated with a trough over Central Asia, which in turn affects the upper-tropospheric circulation over South Asia.
A teleconnection between the North Atlantic and the tropical Pacific is present but rather weak. When using varimax rotation in our M-SSA analysis to reduce spurious correlations in the oscillatory modes, only the 5.9-yr oscillations support a link between the two regions. The energy ratio in case of the 7.8-and 3.6-yr modes, on the other hand, suggests that they are restricted to the North Atlantic and tropical Pacific, respectively.
The teleconnection mechanism between these two regions can be plausibly inferred from previous studies. Chiang and Vimont (2004) have shown evidence that the extratropical atmosphere can influence the tropical Pacific through a meridional mode, analogous to the one in the tropical Atlantic, in which variations in the trade winds in the northern subtropics influence tropical SSTs. Thompson and Lorenz (2004) have identified relationships between the extratropical annular modes and circulation in the tropical troposphere.
The teleconnection between India and the tropical Pacific is supported by the fact that the 2.3-and 2.8-yr oscillations are equally active in both regions and strongly coupled between the two, while a 3.6-yr mode is dominated by the SOI and has but little energy in the Indian core rainfall. The interannual oscillation found in Indian summer rainfall and its teleconnection with ENSO and with tropical Pacific SSTs were the subject of many studies (Sikka 1980;Rasmusson and Carpenter 1983;Parthasarathy et al. 1988;Philander 1990;Kumar et al. 1999;Terray and Dominiak 2005;Sarachik and Cane 2010). In agreement with the recent studies of Penland (1996) and Rajeevan and Pai (2007), we found that the South Asian monsoon is not slaved to the forcing over the equatorial Pacific, but is rather part of a strongly interactive system. Our time-lag analysis pinpointed this to be the case in particular for the 2.3-and 2.8-yr oscillatory modes.
Overall, the approach of synchronized oscillators (Feliks et al. 2010), combined with varimax-rotated M-SSA (Groth and Ghil 2011), appears to be a powerful tool in identifying teleconnections between regional climate modes and helping identify the mechanisms that operate in various frequency bands. This approach should be applicable also to the problems of air-sea interaction and to ocean modes of variability. of the R egion Ile-de-France, while AWR was supported by DOE's EaSM Grant DE-SC0006616 and by ONR's MURI Grant N00014-12-1-0911.

Mathematical and Statistical Details
a. M-SSA Let x 5 fx d (n): d 5 1, . . . , D, n 5 1, . . . , Ng be a multivariate time series with D channels of length N. Each channel is embedded into an M-dimensional phase space by using lagged copies X d (n) 5 [x d (n), . . . , x d (n 1 M 2 1)], n 5 1, . . . , N 2 M 1 1. Note that the resulting matrix X d is constant along the skew diagonals. From this we form the full augmented trajectory matrix X 5 (X 1 , X 2 , . . . , X D ), which has DM columns of length N 2 M 1 1.
The multichannel SSA (M-SSA) algorithm then computes the covariance matrix C 5 X 0 X/N of X and its eigendecomposition; here (Á) 0 indicates transposition. Because of finite-size effects, the sample C may be biased, but effective and accurate estimation methods appear in Ghil et al. (2002b).
Next, one diagonalizes the appropriately estimated covariance matrix to yield a diagonal matrix L that contains the real eigenvalues l k of C and a matrix E whose columns are the associated eigenvectors e k . The e k 's form a new orthogonal basis in the embedding space of X, and the corresponding l k 's give the variance in the direction of e k . The eigenvectors are composed of D consecutive segments e dk of length M: each of which is associated with a channel in X d .
Projecting the time series X onto the eigenvectors, yields the corresponding principal components (PCs) as columns of A.
The pairwise uncorrelated PCs a k are linear combinations of all channels. A way to reconstruct individual channelwise components of the system's behavior in a least squares optimal sense is given by the inverse transformation of Eq. (A2) with a subset of PCs and eigenvectors. Since we are rather interested in the reconstruction of the given channels x d instead of their embedding X d , we average along the skew diagonals of a k e dk 0 , This yields the reconstructed components (RCs; Vautard et al. 1992;Plaut and Vautard 1994). The values of the normalization factor M n and the summation bounds L n and U n for the central part of the time series, M # n # N 2 M 1 1, are simply (M n , L n , U n ) 5 (M, 1, M); for either end they are given in Ghil et al. (2002b). In particular, the time series can be completely reconstructed by the sum of all its RCs, x d (n) 5 å DM k51 r dk (n), as a consequence of å DM k51 a k e dk 0 [ X d . To obtain information about the energy ratio between channels d and d 0 in a given reconstruction r dK (n) 5 å k2K r dk (n), we evaluate here the variance ratio as proposed by Feliks et al. (2010).

b. Varimax rotation of ST-EOFs
The idea of a so-called simple-structure rotation is to find a posterior eigenvector rotation that simplifies the interpretation of the eigenvectors and that reduces mixture effects. There are several ways to quantify the simplicity of an eigenvector's structure (Richman 1986). Varimax rotation attempts to find an orthogonal rotation E* 5 ET that maximizes the variance of the squared elements (Kaiser 1958).
In PCA with M 5 1, the cost functional V 1 being maximized is given by where S is the number of rotated e k *'s, and h * 2 is the corresponding normalization. Kaiser (1958) has given an explicit equation for the sequential rotation of pairs of eigenvectors that shows the algorithm's simplicity of implementation over more sophisticated optimization procedures.
Since the criterion V 1 maximizes the variance over all coordinates, a direct application to ST-EOFs would not only achieve the desirable effect of increasing the spatial variance between the channels, but also increase the temporal variance between the data-adaptive, sine-cosine-like eigenvectors. The latter feature of standard varimax rotation can lead to an undesirable loss of correctly captured oscillatory pairs (Plaut and Vautard 1994).
Recently, Groth and Ghil (2011) have proposed a modification of varimax that maximizes only the variance of the spatial components. Prior to the calculation of the variance, these authors first sum over the temporal part e dk * 2 5 å M m51 e * 2 dk (m) , and so the cost functional V M being maximized becomes with the normalization h * 2 d 5 å S k51 e dk * 2 . In this way, the criterion V M as well as V 1 attempt to bring the spatial component of the S(T)-EOFs either close to one or to zero. In the present paper, we have adopted the modified varimax rotation that takes Eq. (A6) into account. Still, this modified version retains the algorithm's simplicity of implementation: namely, the sequential rotation of pairs of eigenvectors (Groth and Ghil 2011).

c. Annual subsampling
Since we wish to concentrate in the present paper on interannual time scales, we removed all subannual oscillations with a low-pass filter prior to our SSA analysis.
The choice of the filter is critical in order to avoid aliasing effects in the annual subsampling of the monthly time series. An inappropriately chosen low-pass filter cannot only amplitude alter the amplitude of the oscillations but even lead to spurious interannual oscillations in the annually subsampled time series.
To get rid of the strong seasonal cycle in climate records, typical choices include using the annual mean, as well as the seasonal winter or summer means. We compare these customary approaches to low-pass filtering with a Chebyshev low-pass filter. The favored Chebyshev type I filter is very close to an ideal low-pass filter, with practically no distortions of either phase or amplitude in the passband and a highly effective damping in the stopband (Daniels 1974). Figure A1 shows the frequency response of this low-pass filter, together with that of a 12-month as well as a 3-month moving average.
The Chebyshev filter comes quite close indeed to an ideal filter, with a very steep decline at the desired cutoff frequency and a nearly constant magnitude in the passband. The small ripples in the passband can be ignored in comparison to the strong frequency-dependent magnitude response of the two other filters. Although the 12-month filter is more effective in removing high frequencies then the 3-month filter, both filters are strongly frequency dependent and therefore influence any subsequent spectral analysis. In particular, when changing from a monthly to an annual sampling rate, all frequencies above the Nyquist frequency of 0.5 cycles per year are reflected into lower frequencies and the new time series of annually as well as seasonally averaged values are susceptible therewith to spurious oscillations and other aliasing effects.
To demonstrate the different approaches to annual subsampling, we consider the annual mean, along with the summer and winter mean, as well as the July value after Chebyshev low-pass filtering of the given monthly NAO index. In Fig. A2, we estimate the PSD of these annually subsampled records (thick solid) and compare each of them with the PSD estimate of the raw NAO index on the monthly scale (thin solid).
Since the frequency response of the Chebyshev filter has a nearly rectangular shape, we observe practically no difference between the PSD estimates of the monthly and annual time series in the passband (Fig. A2a). The annual mean, on the other hand, matches the monthly PSD only for low frequencies and diverges with increasing frequency (Fig. A2b). Both seasonal-mean values, finally, show large discrepancies over the whole spectral range, with many spurious peaks in the passband (Figs. A2c,d), because of the very weak damping in the stopband. The aliasing in the passband may depend, moreover, on the FIG. A1. Frequency response function of three different low-pass filters: Chebyshev (solid), 12-month moving average (dashed), and 3-month moving average (dotted). The Chebyshev type I filter of order 8 has been set to have its cutoff frequency at 0.45 cycles per year (Daniels 1974). chosen season, with possibly opposite effects in the winter and summer mean. The 7.8-yr oscillation, for example, is completely damped in the summer case and strongly amplified in the winter. Moreover, in the quasi-biennial band we observe even various frequency shifts.
The reconstructed 7.8-yr oscillatory mode is statistically significant in the Chebyshev, annual-mean, and wintermean cases but not in the summer mean. To show the effect of the prefiltering on the SSA analysis, we examine in Fig. A3 the reconstruction of the 7.8-yr mode in each of the annually sampled time series in which it is significant.
In the Chebyshev (Fig. A2a) and annual-mean (Fig. A2b) cases, the 7.8-yr mode is quite similar in PSD amplitude, although a small phase shift emerges (not shown). Strong discrepancies in both amplitude and phase arise between the winter mean (Fig. A2d) and the Chebyshev case (Fig. A2a); these discrepancies can be attributed to the aliasing problem found in the winter mean.
Given the obvious problems with using 3-or 12-month means, we have used in this paper the Chebyshev type I filter for low-pass filtering. Since the results in this case also appear to be independent of the month used for subsampling, the month of January was used throughout, for simplicity's sake.